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In paper a new definition of reduced Pade approximant and algoritlim for its computing is 
proposed. Our approacli is based on tlie investigation of the kernel structure of the Toeplitz 
matrix. It is shown that the reduced Pade approximant always has nice properties which 
classical Pade approximant possesses only in the normal case. The new algorithm allows 
us to avoid Froissart doublets appearance induced by computer roundoff in the non-normal 
O ; Pade table. 

^ ! Pade approximant, Toeplitz matrix, Pade - Laplace method, Froissart doublets. 

Q 

Pade approximants are locally the best rational approximations that can easily be constructed 
^ ! from the coefficients of a given power series. They are closely related to continued fractions, 

^ I orthogonal polynomials and Gaussian quadrature methods [1] . They have been widely used in 

various problems of mathematics, physics, and engineering due to their property to effectively 
solve the problem of analytic continuation of the series beyond its disc of convergence [2]. 

There are many methods available to compute Pade approximants [3]. Some of them 
are implemented in computer algebra systems such as Maple and Mathematica and their 
built-in utilities are frequently used in applied problems. For calculation with floating point 
K*" I numbers in Maple an algorithm due to Geddes [4] based on fraction-free symmetric Gaussian 

elimination is used. Recursive algorithms for computing Pade approximants are also widely 
■ distributed. At first, they allowed to find Pade approximants in the case of a normal Pade 

, table, but later some of them were generalized to the non- normal case; see, for example, [5]. 

^ I It is well known that the Pade table of a rational function R{z) is always non-normal 

since it contains an infinite singular block which elements are identical to R{z). Actually 
the entries inside that block usually differ from the rational function R{z) through the ap- 
pearance of supplementary common roots in the numerators and the denominators, but after 
^ . reducing the common factors, we get the rational function R{z). In any practical calculations 

(because of the computer roundoff and the noise presence in the input data from which Pade 
approximants are constructed), the paired roots in the numerators and the denominators will 
not be rigorously equal. This phenomenon of <pairing» of such zeros and poles got the name 
of Froissart phenomenon and the pairs are known as Froissart doublets [6]. For example, 
their appearance is inevitable in signal processing by using Pade - Laplace method [7]. In 
order to identify doublets, several Pade approximants (besides the desired one) are usually 
calculated. It requires additional coefficients of the Taylor series. It is undesirable, since the 
coefficients are often computed numerically and Pade approximants are known to be very 
sensitive to errors in the coefficients. 

The main purpose of the present paper is to propose a new algorithm for computing Pade 
approximants. This algorithm finds the denominator with minimal degree among all the 
denominators of the Pade approximant and is based on the results concerning the kernel 
structure of Toeplitz matrices [8] . In singular block case it allows to avoid the appearance of 
Froissart doublets induced by computer roundoff. 

The paper is organized as follows. In preliminary Section 2 some definitions and examples 
are given. These examples show that in order to avoid Froissart doublets appearance we 
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should study the kernel of the Toeplitz matrix which gives us denominators of the Pade 
approximant. It is done in Section 3 which contains the main result on the parametrization 
of the set of all the denominators of the Pade approximant. This result allows us to give 
the definition of the modified Pade approximant in Section 4 and establish its properties. 
In Section 5 we propose our algorithm for Pade approximant computing and provide some 
examples in order to show that it allows us to avoid Froissart doublets induced by computer 
roundoff. 



2 Preliminaries 



This section contains some definitions and numerical experiments obtained in Maple and 
Mathematica. We have chosen these packages because they are the most widespread symbolic 
computation systems used in research and applications. 

The examples illustrate how Froissart doublets appear or do not appear when we deal 
with the singular block of a Pade table for a rational function. Note that our new algorithm 
for computing Pade approximants will be proposed in Section 5. Wc will implement it in 
Maple system and solve some of the examples again in order to show that the new method 
allows us to avoid Froissart doublets induced by computer roundoff. 

We start with the classical definition of Pade approximants. 

Definition 2.1 (Pade - Frobenius) 

oo 

Let f{z) be a (formal) power series f{z) — ^ Ckz'^, Ck G C. The (m, n) Pade approximant 

k=0 

corresponding to f{z) is the rational function fm,n{z) — , where Pm,n{z) O'l^d Qm,n{z) 

are polynomials in z such that: 

1- Qm,n{z) ^ 0, diegQm,n{z) < Tl, degPm,n{z) < m, 

2. f{z)Qm,n{^) - Pm,n{^) = r^+n+l^"+"+^ + r^+n+2^™+"+' + • • • • 

In a similar way, we can define the Pade approximant at the point z — a ior the series 

oo 

Ck{z — a)'^. Throughout the paper, we will use Pade approximants at the point z — 0. 

k=0 

Obviously, the coefficients qo, ... ,qn of the denominator Qm,n{z) can be obtained by solv- 
ing the next homogeneous system of linear equations with a Toeplitz n x (n + 1) matrix: 
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Here and further we assume = if A; < 0. From the Condition 2 it follows that, with the 
Qo, . . . , Qn available, the coefficients po, . . . ,Pm of the numerator Pm,n{z) can be obtained from 
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Thus, the {m,n) Fade approximant fm,n{z) always exists. Generally the polynomial 
Qm,n{z) (and hence Pm,n{z)) is found not uniquely, since the rank of the matrix of the 
system (1) can be less than n. Nevertheless, it is easy to show that the rational function 
fm,n{z) is unique. Usually these rational functions are arranged in a double-entry table known 
as the Fade table {/m,n(-2)}m „=o i corresponding to ,f{z). The Fade approximant ,frn,n{,z) 
occupies the nth line and the mth column of the table. It is well known that identical Fade 
approximants can occur only in square blocks of the table. If the Fade table does not contain 
such blocks, it is said to be normal; otherwise it is called non-normal. 

Frovidcd that there is the denominator Qm,n{z) such that (5m,n(0) 7^ 0, we can rewrite 
the Condition 2 as follows 

./X Prn,n{z) _ m+n+l , p m+n+2 , /q\ 

— — j-^ — Itm+n+lZ -\- nm+n+2Z +.... [6) 

In this case, the rational function ^'"'"1^^ is called a Fade - Baker approximant. Note 
that the Fade - Baker approximant does not always exist. 

Let us now come to the second and the main part of the section where we consider 
how Fade approximants can be obtained in such popular mathematical packages as Maple, 
Mathematica and Matlab. 

The built-in utility pade in Matlab approximates time delays by rational linear-time in- 
variant models, i.e. actually it finds only diagonal Fade approximants for an exponential 
function by available for this case explicit formulas [11]. Maple and Mathematica compute 
Fade approximants in the general case and in the examples below we will use their functions 
pade(f, z = a, [m,n]) (in Maple) and Fade Approximant [expr, z, a, m,n] (in Mathematica) 
which compute the {m,n) Fade approximant of the function f{z) at the point z = a. 

To compare the results with the approximated rational functions we will find zeros and 
poles of approximants. 

Example 2.1 Let us find by Maple and Mathematica the diagonal (2,2) Fade approximant 

/2,2 for f{z) = (^^^2i)lz-i) point z = via the following commands. 

In Maple we have: 

_ (z+iyiz-2) . 

J ■ (2+2.1)-(^-l) ■ 

with(numapprox) : p :— pade{f, z — 0, [2, 2]); 

0.9523809523+0.47619047642-0.47619047592-; 
0.9999999999-0.5238095237^-0.47619047592^ 



fsolve(numer(p) , z, complex); 
fsolve(denom(p), z, complex); 



-1.000000000,2.000000001 



-2.100000001, 1.000000000 



In Mathematica we have: 

In[l]:= f := (z + l)*(z - 2)/((z + 2.1)*(z - 1)); 
Inf2j:^ P .- PadeApproximantff {z, 0, {2, 2}}]; F 

niiff9]-= 0.952381+0.476192-0.476192^ 
uaui^J. 1.000000000000000-0.523812-0.4761922 

In[3]:= Roots [Numerator [F] =~ 0, z] 
Out[3]:= z == -l.\\z == 2. 
In[4]:= Roots[Denominator[F] == 0, z] 
Out[4]:^ z -2.1||z 1. 
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As it can be seen, the obtained approximant /2,2 is identical to the function f{z) since 
they have the same zeros and poles. Note that in this case the kernel dimension of the matrix 
of system (1) is equal to 1, i.e. denominator of the Pade approximant is unique. 

It can be verified that /a^s is again identical to the function f{z) (both in Maple and in 
Mathematica), though the kernel is multidimensional. 

However, approximant /4,4 does not coincide with f{z) because of an appearance of sup- 
plementary roots of the denominator and the numerator. 

Example 2.2 Let us find the (4, 4) Pade approximant for the function f{z) — ^^^2 i)lz-i) 
the point z = 0. 

In Maple we have: 

p pade(f z ^ 0, [4, 4j); 

0.9523809524+0.9750566894Z-0 2267573697^^ -0.2494331066^^ 
l.-0.7505668935z^-0.2494331066za 

fsolve(numer(p) , z, complex); 

- 1.909090909, -.9999999999,2.000000000 

fsolve(denom(p) , z, complex); 

-2.100000001, -1.909090908, 1.000000000 

In Mathematica we have: 

In[l]:= P := PadeApproximantff, {z, 0,{4j 4}}]'} P 

niiff1]-= 0-952381+1.1185l2-0.155032z^-0.321159z^-3.33067xl0-i''^^ 
/ V- 1.000000000000000+0.1506242-0.8294652^-0.3211592^+0.24 

In[2]:= Roots [Numerator [P] == 0, z] 

Out[2]:= z == -9.64247 x lQ^^\\z == -1.48273||2 == -l.\\z == 2. 
In[3]:— Roots[Denominator[P] == 0, z] 
Out[3]:^ z = -2.1||z -1.482731 1^ 1. 

Thus, in Example 2.2 Froissart doublets induced by computer roundoff appear. Note 
that doublets obtained in Maple and Mathematica are different (they are in bold) . The root 
—9.64247 X 10^^ of the numerator (obtained in Mathematica) is caused by its small leading 
coefficient -3.33067 x 10^^^ 

Note that for /a s and /4 4 the denominator is not unique and, as we have seen. Maple 
and Mathematica may choose not the best one. 

Let us take another rational function and find its approximation in its infinite singular 
block. As it can be seen below. Maple does not produce Froissart doublets, but they appear 
in Mathematica. 



Example 2.3 Let us find the (2, 3) Pade approximant for the function f{z) = (^^_^2){l-2 01) 
at the point z = 0. 
In Maple we have: 



z+l.Ol 



■ (z+2)-(2-2.01) ■ 

with(numapprox): p :— pade(f z, [2, 3]); 



-0.2512437812-0.24875621902 



1.000000000+0.0024875619922-0.2487562188x2 

fsolve(numer(p) , z, complex); 

-1.010000000 

fsolve(denom(p), z, complex); 

-2.000000001, 2.010000000 
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In Mathematica we have: 

In[l]:=f:= (z + 1.01)/((z + 2)*(z - 2.01)); 

In[2]:^ P PadeApproximantff, {z, 0, {2, 3}}J; P 

niif[9h= -0.251244-0.1220692+0.125433^^ 

/^V- 1.000000000000000-0.501752^-0.2500112^+0.1254332^ 

In[3]:= Roots [Numerator [P] == 0, z] 
Out[3]:= z == -1.01||^ == 1.98318 
In[4]:— Roots [Denominator[P] 0, zj 
Out[4]:^ z -2. 11^ 1.983181 1^ 2.01 

The following conclusions may be drawn from the examples. When the kernel dimension 
of the matrix of system (1) is equal to 1 (as in Example 2.1), we have one denominator 
of the Pade approximant and Proissart doublets can not appear. But they can sometimes 
appear when the kernel dimension is greater than 1. If an improper denominator was chosen, 
then supplementary (also called artificial) roots appear. These roots may have corresponding 
pairs among the roots of the numerator and the cancelation of the common factors reduces 
the Pade approximant to the reduced form. But in the presence of computer roundoff, the 
paired roots in the numerator and denominator will not be rigorously equal and this can 
cause the appearance of Proissart doublets (as in Examples 2.2, 2.3). Moreover, computer 
roundoff may produce artificial roots with the great absolute value (as in Example 2.2), since 
the vanishing leading coefficients are small but not equal to zero. 

To avoid these problems, we should choose the denominator with the minimal degree. As 
will be shown further, this denominator always exists. In order to find it, we have to study 
the kernel structure of the Toeplitz matrix of the system (1). It will be done in the following 
section. 



3 Parametrization of the denominator set 

In this section we study the structure of the set of all denominators for a (m, n) Pade ap- 
proximant, i.e. the structure of the kernel kerT^+i of the Toeplitz matrix 

Tm+l — i = m+l,...,m + n 

j = 0,l,...,n 

from the system (1). 

Pirst of all, let us prove that the minimal degree denominator exists and establish its 
properties. 

Proposition 3.1 There exists the denominator Q^^ni^) with the minimal degree among all 
the denominators Qm,n{z). This denominator is unique up to a constant factor. 

Let P!^^n{z) be the numerator of the Pade approximant corresponding to (5^ „(z). Then 
the polynomials -P^,j(-z) and Q^,ni'^) '^^^ have common non-zero roots. 

Proof. Since all the polynomials Qm,n{z) satisfy the condition (legQm,n{z) < n, there exists 
the denominator with the minimal degree d. 

Suppose that there are two denominators with the degree d: 

Ql^^iz) =Bdz'' + ... + Bo, Ql^ni^) - Bdz'' + ... + Bo, B^^O, B^^O. 

Let us introduce Q{z) = BdQ^„^^{z) — BaQ^^i^). This is the denominator of the Pade 
approximant with the degree less than d. It is possible if and only if Q{z) = 0. Thus, Q^ni^) 
and Qm,n(^) linearly dependent. The uniqueness is proved. 



5 



Let P^,ni^) be the numerator corresponding to Q^,^{z). Let us suppose that zq is the 
common non-zero root of P^ni^) Qmni^)'- 

By Definition 2.1 we have: 

Then „ (2;) is the denominator with the degree less than d. We have a contradiction since 
d is the minimal degree. 

■ 

Now we describe the denominator set {Qm,n{z)}, i.e. the kernel structure of T^+i- This 
result will be used to construct the algorithm for obtaining the minimal degree denominator. 

To find any denominator of the Pade approximant, we need the sequence of the Taylor 
coefficients Cm_„+i, . . . , c„i+n of f{z). The sequence consists of the entries of matrix T^+i. Let 
us introduce notions of indices and essential polynomials for Cm-n+i, ■ ■ ■ , c^+n- The notions 
were given in more general situation in the paper [8]. Here we formulate them specially for 
our case. 

It is natural to include T^+i in the family of Toeplitz matrices 



Cfc+l Cfc • • • Cm-n+2 



m — n + l<k<m + n, (4) 



which are constructed from the elements of Cm-n+i, ■ ■ ■ , Cm+n- 

Consider the sequence of the spaces ker T^. It is more convenient to deal not with vectors 
Q — {Qo, Qi, ■ ■ ■ , Qk-m+n-iY £ ker but with their generating polynomials Q{z) = qo + cji z + 
• • • + qk-m+n-i Instead of the spaces kerT^ we will use the isomorphic spaces Afc 

consisting of the generating polynomials. 

To do this, we introduce a linear functional a by the formula: 

a{z^} = C-j, —m — n < j < —m + n — 1. 

In the theory of orthogonal polynomials this functional is called the Sticltjcs functional. 

Denote by Mk {m — n + 1 < k < m + n) the space of polynomials Q{z) with the formal 
degree k — m + n — 1 satisfying the orthogonality conditions: 

a{z~'Q(z)} ^0, k,k + l,...,m + n. (5) 

It is easily seen that N'k is the space of generating polynomials of vectors in kerT^. For 
convenience, we put Mm-n = and denote by Afju+n+i the (2n + l)-dimensional space of all 
polynomials with the formal degree 2n. 

Let dk be the dimension of the space J^k and = dk — dk-i {m — n+ 1 < k < m + n+ 1). 
The following fact is crucial for the further considerations. 

Proposition 3.2 For any non-zero sequence Cm-n+i, ■ ■ ■ , Cm+n the following inequalities 

= Ajn~n+1 < < • • • < ^m+n < ^m+n+l — 2 (6) 

are fulfilled. 
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Proof. It follows from Definition 2.1 that A4 zNk are subspaces of A4+i {m—n < k < m+n) 
and 

Hence, by the Grassman formula, we have: 

dim(7Vfe + z^^k) = 24 - 4_i. (7) 

Let hk+i be the dimension of any complement T-Lk+i of the subspace + zAf^ in the 
whole space A4+i. Prom (7) we have 

hk+i = Afc+i — Ak, (8) 

i.e. Afe+i > Afe. By definition d^^n = and d^-n+i is also equal to zero. Hence, Ajn_n+i — 0. 
In a similar manner we can prove that Am+n+i = 2. ■ 
It follows from the inequalities (6) that there exist integers /Ui < ^2 such that 

^m-n+l — • • • — = 0, 

A^,+i = ... = A^, =1, (9) 

If the second row in these relations is absent, we assume /ii — fj,2- Really in our case /ii < ^2 
as will be shown later. 

Definition 3.1 The integers /J^i, defined in (9) will he called the essential indices (briefly, 
indices) of the sequence Cm-n+i, ■ ■ ■ , Cm+n- 

Proposition 3.3 Let x = rankT^. Then the indices /Ui,//2 are found by the formulas: 

III — m — n + X, fj,2 — TTi + n — 

Proof. It follows from the definition of A^ that ^ A^ = dm+n+i — dm-n — 2n + l. 

k=m—n+l 

On the other hand, from the relations (9) we have 

m+n+l 



Afc = 1 • (/X2 - Ail) + 2 • (m + n + 1 - //a)- 



fe=m— n+l 

Hence, /Ui + /i2 = 2m + 1. It means that /ii < m < ^^2- 

m m 

In a similar manner we obtain ^ A^ = dm = n — x and ^ Ak = m — fii. 

k=m—n+l fe=m— n+l 

Thus, Hi = m — n + x. 

Since Hi + H2 = 2m + 1, we get iJ,2 = rn + n — 

■ 

Now we can describe the structure of the kernels of the matrices Tk- It follows from (8) 
and (9) that hk+i :— dimHk+i iS k — /ij {j — 1,2). In that case 4+i — 1- Therefore, 
for k 7^ /ij 

Nk+i=Nk + zNk, (10) 

and for k = /ij 

Nk+i = {Mk + zMk)+Hk+i. (11) 
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Definition 3.2 Any polynomial Qj{z) that forms a basis for one- dimensional complement 
will he called the essential polynomial of the sequence Cm-n+i, ■ ■ ■ , c^+n corresponding 
to the index /ij, j — 1,2. 

It can be shown (see [8], theorem 4.1) that integers fii, are the indices and polynomials 
Qi{z) e Q2{z) e M^^+i are the essential polynomials iff 

ao := a{2-™-"-i[Q2(0)Qi(2) - QmQ2{z)]} ^ 0. (12) 

In the following theorem the structure of the kernels of matrices (m— n+1 < k < m+n) 
is described in terms of the indices and the essential polynomials. 

Theorem 3.1 Let /Ui,A<2 be the indices and let Qi{z),Q2{z) be the essential polynomials of 
the sequence c^-n+i^ ■ ■ ■ ■> Cm+n- 
Then 

{0, m — n-\-l<k< III, 

{qi{z)Qi{z)}, Hi + l<k<fX2, 

{qi{z)Qi{z) + q2{z)Q2{z)}, H2 + i<k<m + n, 

where qi{z), q2iz) are arbitrary polynomials of the formal degree k — /ij — 1. 

Proof. From (9) we have dk = for k E [m — n + 1, fii), i.e. A4 = 0. 

Let fjii + I < k < ijL2- It follows from (10) and (11) that the polynomials 

{Ri{z), zR,{z), z'-^^-'R.iz)} (13) 

generate the space Afc. The number of these polynomials is equal to k — /ii. 

k 

Prom the definition of and the relations (9) we have dk— ^ Aj — k — 

j=m—n+l 

Thus, the number of the polynomials in (13) is equal to the dimension of the space A^- 
Therefore, these polynomials form a basis of J\fk. 

The case ii2 + i'^k<m-\-n can be considered in a similar manner. 

■ 

Now we apply the previous theorem to the case k = m-\- 1 and obtain the main result of 
this section on the parametrization of the denominator set. 

Theorem 3.2 Let Qm,n{z) be an arbitrary denominator of the {m,n) Fade approximant for 
00 

the series f{z) = ^ Ckz''. Form the sequence {cm-n+i, ■ ■ ■ ,Cm+n} which is necessary to find 

k=0 

Qm,n(^Z^ . 

Let III be the first index, let Qi{z) be the first essential polynomial of this sequence. 
Then the denominator set is 

{Qm,n{^)} = {?i(^)<5i(^)}> 

where qi{z) is an arbitrary polynomial with the formal degree m — fj,i. Thus, Qi{z) is the 
denominator Q^ ,^{z) with the minimal degree. 

Remark 3.1 It follows from Qm,n{,z) = q{z)Qi{z) that the denominator Qm,niz) such that 

Qm,n{Q) 7^ (Baker's condition) exists iff Qi{0) ^ 0. 

Hence, the Fade approximant with the minimal degree denominator is the Fade - Baker 
approximant if the latter exists. 

In the following section we use this parametrization in order to modify the definition of 
the Fade approximants. 
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4 Reduced Pade approximant and its properties 



Definition 2.1 ignores the non- uniqueness (in general case) of a denominator of a Pade approx- 
imant. It does not matter mucii for tlie normal Pade table but in the singular case common 
factors of the denominator and the nominator can not be cancelled because of roundoff errors. 
This is one of reasons why Proissart doublets can appear. 

The main aim of this section is to modify the classical Definition 2.1. We would hke a 
Pade approximant to exist always, to have a unique denominator (up to a constant factor) 
and to coincide with the Pade - Baker approximant if the latter exists. 

Let us add the minimality requirement to Definition 2.1. 

oo 

Definition 4.1 Let f{z) be a (formal) power series f{z) = ^ CkZ^, Ck G C. The {m,n) 

k=0 

Pade approximant corresponding to f{z) is the rational function fm,n{z) — , where 

Pm,n{z) and Qm,n{z) are polynomials in z such that: 

1- Qm,n{z) ^ 0, (\egQm,n{z) < U, dcg Pm,n{z) < TU, 

2. f{z)Qjn,n{z) — Pm,n{z) — rm+n+lZ"^~^^~^^ + rm+n+2Z"^~^"'~^'^ + . . . . 

3. The polynomial Qm,niz) has the minimal degree among all polynomials satisfying 1, 2. 

Theorem 3.2 on parametrization of the denominator set gives the constructive method 
for finding the denominator with the minimal degree. 

The next theorem shows that the Pade approximant has the desired properties. 

Theorem 4.1 For any power series f{z) the following statements are fulfilled. 

1. Any {m,n) Pade approximant exists and is unique. 

2. The denominator Qm,n{z) of the Pade approximant is unique up to a constant factor 
and Qm,n{z) is the first essential polynomial Qi{z) of the sequence {cm-n+i, ■ ■ ■ , Cm+n}- 

3- Pm,n{z) and Qm,n{z) have not common non- zero roots. 

4. The Pade - Baker approximant exists iff Qi{0) 7^ 0. The Pade approximant from 
Definition 4-i is the Pade - Baker approximant if the latter exists. 

5. If Qi{z) has the root z — of order 5m,n > 0, then 

f{z) - fm,n{^) = + + . . . , A^O. 

(^m,n 'is called the deficiency index of Pade approximant.) 

Proof. Statements 1-4 evidently follow from Proposition 3.1, Theorem 3.2 and Remark 3.1. 

Prove the last statement of the theorem. Let z = be the root of order 5m,n of Qi{z). 
Since Qi(0) = 0, the number uq from formula (12) is 

<To = -a{z-^-^-^Q,{z)}Q2{Q)^Q. 

Here Q2{z) is the second essential polynomial. 
Thus ct{z-'^-"-iQi(z)} 7^ 0. 
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It is easy to see that the number a{z " ^Qi{z)} is the coefficient at z^^'^^^ in the 
power series f\z)Qi{z). As it is not zero, we have 



f{z)Q,{z) - P,{z) = + . . . , A = a{z-^-^-^Q^{z)] + 0. 

Here P\[z) is the numerator corresponding to the denominator Q\{z). 
After dividing this equation by Q\{z)i we get 

where A ^ 0. 

It means that the deficiency index of the Pade approximant fra,n{^) coincide with the 
multiplicity Sm,n of the root 2; = of the first essential polynomial Qi{z). ■ 

Since the numerator and the denominator of Pade approximant from Definition 4.1 have 
not common non-zero roots, we will call it reduced Pade approximant. Note that the numer- 
ator and the denominator of the reduced Pade approximant can not be coprime because of 
their common zero roots. 

Remark 4.1 It is impossible to improve Definition 4-i in such a way that the denominator 

and the numerator are always coprime. 

It is easy to see that if Qm,n{^) = z^Qm,ni^), then Pm,n{z) = z^P^„ „{'^). Hence /m,„(z) = 

"I''" . , . However the polynomials P^^{z), Q^^[z) can not be considered as the numerator 

and the denominator of fm,n{^)} ^'ince the Frobenius Condition 2 does not fulfilled: 

f{z)QlA^) - P°,„(^) = + . . . . 

To determine 5 (the deficiency index of fm,n{z))i should know how many low order 
coefficients of Qi{z) are equal to zero. Note that vanishing coefficients of the denominator 
and the numerator can be non-zero because of roundoff errors. In particulary, it can cause 
the appearance of roots of the denominator and/or the numerator with the great absolute 
value (as in Example 2.2). 

Thus, we face the problem of finding vanishing coefficients of the denominator and the 
numerator. The results of Theorem 4.1 allow to solve this problem. This will be done in the 
following two theorems. 

Recall that x = rankT^, /ii = m — n -|- x, matrix T^^+i has the size (2n — x) x (x-|- 1), 
and rankT^j+i = x. The vector which is a basis for one-dimensional space kerT^^+i gives 
coefficients of the first essential polynomial Qi{z) which is the minimal degree denominator 
of the Pade approximant. 

Further we will denote by T^^+i the matrix which is obtained from the matrix T^^+i by 
deleting of the fcth column, k — l,...,x-|-l. 

Theorem 4.2 Let Qi{z) = qQ + qiz + - ■ ■ + q^z'^ be the denominator with the minimal degree. 
Then qk = ^ iff rank T^^^ = x — 1, = 0, . . . , x. 

Proof. Let qu = 0. We have YBnkTl^lXi < Suppose rankT^^:|:i^ = x. Then T^^+i^ has the 
trivial kernel, i.e. go = 9i = • • • = Qk-i = Qk+i — ■ ■ ■ — q^ — 0- Thus, Qi{z) = 0. Since it is 
impossible, we have rankT^J^^J^ < x — 1. 

Assume that rank T^^^ J'' < x — 1. Then T^J^+i^ has a multidimensional kernel, which can 
be embedded in the kernel of T^^+i by a natural way. But ker T^^+i is one-dimensional, hence 
our assumption is false and rankT^J^^^'' = x — 1. 
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On the other hand, if rank 7"('^'+^) 



X — 1, then kerT^f^il^ is one-dimensional and after 



>i+i — ^"^^ -^//i+i 

the natural embedding in one-dimensional space ker T^^+i we get qk — 0. 

■ 

Denote by T^^+i the matrix which is obtained by inserting the row (cfe Cjfc_i . . . Cjt_^) at 
the beginning of the matrix T^j+i, k = 0, . . . ,m. 



Theorem 4.3 Let Pi{z) — po + Piz -|- • • • -|- PmZ^ be the numerator corresponding to the 
denominator Qi{z). Then Pk — rank T^^'_,_^ = k — 0, . . . ,m. 

Proof. It follows from (2) that Pk = CkQa + c^-iqi + ■ — h Ck-^Qx, k = 0, . . . ,m. Hence, 



T. 



[k] 







If pk = 0, then matrix Tjj'^_^_i has a nontrivial kernel, hence rank r^^'_,_j^ < x + 1. Since 
rankT^^+i = x, we have TankTjj'^_^_i = x. 

On the other hand, if rankT^^^'^^ = x then the inserted row (c^ Ck-i ■ ■ ■ Ck-x) is a linear 
combination of the rows of T^i+i- Therefore, CkQa + Ck-iqi + • — h Ck-^Qx — 0, i.e. Pk — 0. 

■ 

From Theorems 4.2, 4.3 we can obtain the multiplicities of 2; = 0, 2; = oo as the roots of 
Qi{z), Pi{z), i.e. the deficiency index of fm^niz) and the degrees of Qi{z), Pi{z). 
In particularly, the next result on the deficiency index is now evident. 

CoroUciry 4.1 The number S is the deficiency index of fm,n{z) 'iff 

rankT^J^i = . . . = rankT^^^^^ = x - 1, rankT^f+i = x. 



We would like to end this section with the following remark. As is well known in the case 



of norma/ Pade approximant /m,„(2;) 



n,,n (-2) 



the degrees of Pm,n{z) and Qm,n{z) are equal 



to m and n, respectively, and 6 = 0. Thus, in the normal case degrees of the numerator and 
the denominator and the deficiency index 6 are known. Our modified Definition 4.1 and the 
results of this section allow us to determine them not only in the normal case, but also in 
the non-normal one. 



5 Algorithm 

In this section we present our algorithm for computing a reduced Fade approximant with the 
minimal degree denominator. As we have seen in the previous section, in order to realize the 
algorithm, we have to find the rank of matrix T^, and the null space of T^i+i. Moreover, 
in order to delete vanishing coefficients of the denominator and the numerator we have to 
find ranks of matrices TJ^^+i, k = 1, . . . , x + 1, and TJ^^+i, k = 0, . . . ,m. In practice, due to 
rounding and measuring errors and finite computer precision, the elements of these matrices 
are perturbed with error, so that we have to determine the rank and the null space of the 
original matrix from the perturbed matrix. 

For computer calculations we need to define the numerical rank and the numerical null 
space of A. Recall these concepts which are crucial for the algorithm. 
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The definition of tlie numerical rank was first given by Golub, Klema and Stewart [9]. We 
will use the simplified definition (as in [10], [11], p. 72). The numerical e rank oi an M x N 
matrix A with respect to the threshold £ > is defined as the smallest rank of all matrices 
within a 2-norm distance e of A. Namely, 

rank [A, e) = min{rank B : \ \A — B\\2 < s} . 

The numerical e rank may be characterized in terms of the singular value decomposition 
(SVD). According to Stewart [12], «;the singular value decomposition is the creme de la creme 
among rank-revealing decomposition». 

Let us recall the definition of the SVD. Every MxN complex matrix A can be represented 
in the form 

A = UEV", 

where U, V are unitary matrices means the Hermite conjugation), is an M x N matrix 
in which the upper N x N block is a diagonal matrix with all entries real and sorted in 
descending order. The diagonal entries ui > . . . > (Tjv of the matrix E are called the singular 
values. 

In the case of exact calculations, we have o"i > . . . > ct^ > cr^+i = . . . = = 0, where r 
is the rank of A. In practice, we have 

<Ji > . . . > ar > S > ar+l > C"jv > 0. 

Let us denote Ar — LfErV^ with = diag{(Ti, . . . , (7^, 0, . . . , 0}; then 1 1 A — 1 12 = c^+i and 
rank (A, e) = rank A,. = r [9]. Moreover, A^ is the nearest matrix to A (with respect to the 
2-norm) with rank r. Therefore the null space of A^ is called the numerical null space of A 
within e. The null space of Ar is spanned by {vk+i, ■ ■ ■ , v^}, where Vk is the kth column of 
the matrix V. The ratio 7 = is called the numerical rank gap. The numerical rank of 
the matrix A may be estimated reliably in the case of «:well defined numerical rank» that is 
when A has a single well-determined gap between large and small singular values. 

Thus, in terms of singular value decomposition, numerical e rank is the number of singular 
values greater than the given threshold e. There is no uniform threshold for all applications. 
The user must make a decision on the threshold e based on the nature of the applications. 

The rank(A,tol) function in Matlab returns the number of singular values of A that are 
larger than the threshold (tolerance) tol. The default tolerance is max(M, A'")*eps(norm(A)). 
Here norm(A) indicates the Euclidean norm of A and eps(norm(A)) is approximately 2.2 x 
10~^^ times norm(A). This choice is usually a good choice if the errors in the matrix elements 
are due to computer arithmetic and if there is a sufficiently large gap in the singular values 
around this tolerance. 

Although the SVD is the most widely used method for determination of the numerical 
rank and the numerical null space, there are alternative methods like URV decomposition, 
LU decomposition or QR factorization with column pivoting. 

Now we can present our algorithm for computing (m, n) reduced Fade approximant 
fm,n{z) for f{^) at the point z — a. 
Algorithm. 
Initialization: 

f{z) := the approximated function; 

(m, n) := the order of the Fade approximant; 

a :— the center of expansion; 

(i := 0, if we want to delete vanishing coefficients of numerator and the denominator, 
else d :— 1. 



12 



1. Compute the Taylor coefficients cq, . . . , Cm+n of f{z) at the point z = a. 

2. Form the Toephtz matrix = ||cj_j+m|| i = i,...,n + i, 

j = 1, . . . ,n. 

3. Determine its rank x = rankT^ and the index /ii = m — n + x. 

4. Form matrix T^^+i = i = i,... 

J = 1, . . . , n — m + /ii + 1. 

5. Find a basis for its one-dimensional kernel kerT^j+i. 

The obtained vector (^oi 9i) • • • i 9n-m+^ti)^ is a vector formed from the coefficients of the 
minimal degree denominator Qi{z) — ^ qk{z — of the Fade approximant. 

fc=0 

6. li d = 0, then for /c = 0, . . . , x, do: 

a. Form the matrix T^^^^ . 

b. Find its rank: rankT^^^}'*. 

c. If rank r^^';J:J^ = x - 1, then replace in (go, gi, • • • , gn-m+^^J^ by zero. 
Else step 6 should be omitted. 

7. For = 0, . . . ,m do pfc = Ckqo + Ck-iqi-\ \-Ck-xqx and form the vector (po, • • • ,]?m)^ 

consisting of the coefficients of the numerator Pi{z) = ^ pk{z — a)'' corresponding to Qi{z). 

k=0 

8. li d = 0, then for A; = 0, . . . , m, do: 

a. Form the matrix T^j+i. 

Ik] 

b. Find its rank: rankT^^_,_p 

c. If rankT^^'_,_^ = x, then replace pk in (po, . . . ,Pm)^ by zero. 
Else step 8 should be omitted. 

End of Algorithm. 

Output: fm,n{z) = 

Comments: 

1) Step 1 consisting of computing the Taylor coefficients is very important for the algo- 
rithm. It is well known (see, for example, [3]) that accuracy in the given coefficients Cfc is 
essential for Fade approximants. Usually most computing effort goes into calculation of the 
coefficients rather than Fade approximants, and so the coefficients should be calculated with 
the greatest possible accuracy. The way of their calculating depends essentially on the prob- 
lem under consideration. For example, in the Fade - Laplace method coefficients are found 
from experimental data by numerical integration. So the accuracy in the coefficients depends 
not only on the chosen method of integration but also on the quality of the experiment. 

2) On step 3 we find the essential index pLi according to Froposition 3.3. The ffist 
essential polynomial, which is the minimal degree denominator of the Fade approximant, 
is found on step 5 according to Theorem 3.2. Theorems 4.2, 4.3 are used on steps 6, 8, 
respectively. The output approximant fm,n{z) possesses properties listed in Theorem 4.1. 

We have implemented the algorithm in Maple (procedure ReducedPade(f,m,n,a,dJ) and 
now we would like to repeat Examples 2.2, 2.3 from Section 2. 

Example 5.1 Let us find by ReducedPade the diagonal (4,4) Pade approximant fa for 
~ (i+2j)(^-i) pom^ z — via the following commands. 

J ._ (z+l)-(z-2) . 



(z+2.1)-(2-l) ■ 

app :— ReducedPade(f 4, 4) 0, 0); 

0.7773220744+0.3886610375z-0.388661037l2^ 

0.8161881781-0.4275271406546188322-0.388661037357416250z2 

fsolve(numer(app) , z, complex); 
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-0.9999999998, 2.000000001 
fsolve(denom(app), z, complex); 

-2.099999999, 1.000000000 

Example 5.2 Let us find the (2,3) Pade approximant for the function f{z) — (_^_^2)(z-2 oi) 
at the point .^' = 0. 



-1.01 



f ■ (^+2)-(z-2.01) • 

app := ReducedPade(f, 2, 3, 0, 0); 

-0.2438127455-0.2413987580Z 

0.9704230069+0.002413987817669799312-0.2413987580543237872^ 

fsolve(numer(app), z, complex); 

-1.010000000 

fsolve(denom(app), z, complex); 

-1.999999999, 2.010000000 

As it can be seen, the obtained approximants are identical to the approximated functions 
and Froissart doublets do not appear. Note that in the given examples the value of the 
parameter d was equal to zero. Hence, the vanishing coefficients of the numerators and the 
denominators were deleted. The next example shows that such deleting is desirable. 

Example 5.3 Let us find the (3,7) Pade approximant for the function f{z) — ^a^^^2^^ cl^ 
the point z = 0. 

t ._ 2 + 1.01 . 

J ■ 24+3^2-4.01 • 

app :— ReducedPade(f 3, 7, 0, 0); 

0. 1977728838+0. 19581473642 

-0.7852170931+0.58744421453119533022+0.1958147376729382782-* 

fsolve(numer(app), z, complex); 

-1.010000000 

fsolve(denom(app) , z, complex); 

-1.000999098, -2.0004997387, 2.0004997387, 1.000999098 

app :— ReducedPade(f 3, 7, 0, 1); 

(0.1977728838 + 0.1958147377Z- 1.48862691 10-"z3)/ 
(-0.7852170931 - 5.23710724786852211 lO'^^; + 0.587444214531194664^2 _ 1.13198619922094679 lO'V + 
0.1958147376729380560^ - 1.67404646568558580 lO-^o^^) 

fsolve(numer(app), z, complex); 

-1.146906047 10^, -1.009999994, 1.146916147 10^ 

fsolve(denom(app), z, complex); 

-1.000999095, -1.728850777 lO'^ - 2.0004997387, -1.728850777 lO'^ + 2.0004997387, 

1.000999101, 1.169709095 10^ 
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6 Conclusions 



This work was motivated by our intention of avoiding the Froissart doublets appearance in 
signal processing by Fade - Laplace method. To do this, we have modified the classical 
definition of Fade approximant by adding the requirement of the minimal degree of its de- 
nominator. It turned out that this new definition of reduced Pade approximant allows us to 
avoid Froissart doublets induced by computer roundoff. 

The reduced Pade approximant can be easily obtained by the proposed algorithm and 
always has nice properties which classical Fade approximant possesses only in the normal 
case: the denominator is unique up to a constant factor, the numerator and the denominator 
have not common non-zero roots, their degrees and the deficiency index are known exactly. 
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